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1.0  INTRODUCTION 


This  final  rsporc  concludes  our  study  of  large-scale  vectorized  numeri¬ 
cal  analysis  applied  to  time  domain  seismic  wave  phenomena  in  filled  basins. 
Applications  include  calculations  of  waves  from  simple  surface  or  buried 
sources  in  a  variety  of  idealized  2-D  Basin  and  Range  models  (36,000  to 
120,000  nodes)  and  one  large  3-D  model  (400,000  nodes)  from  Yucca  Flat, 

Nevada  Test  Site.  Analysis  is  based  on  an  explicit,  finite  element,  elastic 
wave  solver  designed  for  fully  vectorized  execution  (vector  processing  of 
code  loops)  on  the  CRAY-1.  This  combination  of  an  explicit  solver  and  the 
new  generation  of  supercomputers  permits  full  wave  solutions  in  inhomogeneous 
models  10  to  100  times  faster  and  larger  than  feasible  using  conventional 
solvers  on  scalar  mainframe  computers  (e.g.,  CDC  7600). 

The  study  is  motivated  by  interest  within  the  U.S.  Air  Force  in 
calculations  of  high  frequency  wave  fields  in  realistic  models  of  basins 
and  basin-like  geology.  Well  controlled  models  generally  include  smooth 
wave  speed  variations  within  each  material,  as  well  as  discontinuities, 
curved  interfaces,  pinch-outs,  faults,  sharp  corners,  irregular  topography 
and  surface  layers.  Typical  analytical  methods  require  significant  approxi¬ 
mations  to  make  these  problems  tractable.  For  example,  simplifying  assump¬ 
tions  on  the  wave  field  often  involve  high-frequency  asymptotics  and  ray 
tracing  which  usually  neglect  diffractions  and  surface  waves,  or  sometimes 
the  replacement  of  elastic  by  acoustic  media  in  discrete  numerical  models. 
Physical  idealizations  can  lead  to  plane,  layered  models  which  admit  semi- 
analytic  half space  solutions;  or  piecewise  constant  material  properties  and 
smooth  Interfaces  amenable  to  boundary  Integral  methods;  or  discrete  approxi¬ 
mations  but  with  reduced  model  size  to  fit  into  fast  computer  memory.  However, 
with  the  computational  power  afforded  by  explicit  methods  and  supercomputers, 
these  idealizations  are  often  unnecessary  and  more  complete  or  accurate 
solutions  are  available.  In  any  case,  large-scale  calculations  are  useful 
in  general  modeling  applications  and  numerical  experiments,  as  well  as  the 
evaluation  of  modeling  Idealizations. 

Research  goals  of  the  study  are  two-fold.  First  is  to  perform  large- 
scale  calculations  and  seismic  interpretation  for  a  number  of  interesting 
geologic  models.  Second,  in  support  of  the  calculations,  is  to  study 


theoretically  the  numerical/computational  basis  for  discrete  hyperbolic  wave 
solvers  (explicit  finite  difference/finite  element  methods)  and  to  Increase 
their  utility  in  seismic  applications,  and  their  compatibility  with  super¬ 
computer  architecture.  This  latter  goal  includes  an  evaluation  of  present 
capabilities  on  available  CRAY  hardware. 

A  previous  report  addresses  these  goals  in  the  context  of  2-D  modeling 
of  some  typical  Basin  and  Range  geology.  An  overview  of  results  is  presented 
in  1.1  below.  The  present  report  is  concerned  principally  with  3-D  calcula¬ 
tions  for  a  model  of  the  recent  COALORA  test  in  Yucca  Flat,  Nevada  Test  Site. 
Due  to  the  complexities  associated  with  3-D  modeling  and  interpretation,  the 
Yucca  Flat  site,  suggested  by  Stump  (1983),  was  chosen  over  other  candidate 
sites  because  an  excellent  geophysical  data  base  exists,  incorporating  seismic, 
gravity  and  well  log  data.  The  geophysical  data  base  and  models  are  discussed 
in  Section  2  and  the  2-D  and  3-D  finite  element  calculations  are  described  in 
Section  3.  A  theoretical  analysis  of  the  explicit  wave  solver  algorithm  is 
contained  in  Section  4.  A  brief  discussion  and  conclusions  complete  the  report 
in  Section  5. 

1.1  BACKGROUND 

The  following  is  a  summary  of  results  from  an  interim  report,  Wojcik, 
et  al.  (1982),  containing  a  numerical  analysis  of  2-D  basins  typical  of  the 
Basin  and  Range  province.  Calculations  were  performed  to  study  the  effects 
of  basin  edge  geometry  on  the  reflection,  transmission  and  diffraction  of 
body  and  surface  waves  excited  by  surface  sources.  Normal  and  tangential 
surface  tractions  were  applied  on  the  model  centerlines  and  synthetic  seis¬ 
mograms  generated  by  convolving  a  2  Hz  wavelet  with  numerical  Green's 
functions  recorded  at  closely  spaced  stations  over  the  surface.  Seismic 
velocity  functions  used  piecewise  linear  and  piecewise  constant  approximations 
to  a  suite  of  functions  compiled  by  Battis  (1981)  from  available  Basin  and 
Range  geophysical  data.  Fig.  1.1-1. 

Symmetric  basin  models  with  three  types  of  flank  structure  were  examined 
first.  The  models  and  synthetic  seismograms  for  normal  and  tangential  sur¬ 
face  loading  are  illustrated  in  Fig.  1.1-2  and  1.1-3,  respectively.  These 
and  other  calculations  showed  a  number  of  quantitative  differences  between 
a  steep  faulted  flank,  an  echelon  faulted  flank  and  a  dipping  flank.  Basin 


edges  are  efficient  transmitters  of  Rayleigh  waves,  apparently  independent 
of  edge  details;  however,  geometry  has  a  strong  effect  on  Rayleigh  wave 
reflection.  Shear  waves  outside  the  basin  depend  markedly  on  source  type 
and  edge  geometry  due  to  mode  conversions  there.  The  strong  shear  phase 
generated  by  the  tangential  source  is  converted  primarily  to  S-waves  at  the 
edge,  whereas  the  shear  phase  from  the  normal  source  converts  most  of  its 
energy  to  Rayleigh  waves. 

A  model  of  two  basins  separated  by  a  mountain  was  also  examined.  This 
basin-range-basin  model.  Fig.  1.1-4,  extended  the  above  results  for  a  single 
basin  and  provided  Insight  into  the  type  and  strength  of  signals  communica¬ 
ted  between  basins.  The  finite  element  Idealization  had  120,000  elements, 
required  20  minutes  of  central  processor  time  (overnight  cost,  $100.)  to 
propagate  the  Rayleigh  wave  across  it,  and  demonstrated  that  large-seal' 
full  wave  field  calculations  are  practical  on  the  CRAY-1. 

The  major  result  of  the  finite  element  theoretical  work  was  that, 
the  context  of  Cartesian  (rectangular)  finite  element/finite  difference 
grids,  the  governing  equations  on  the  element  level  can  be  greatly  simplified 
by  linear  transformation  and  formulated  so  as  to  admit  selective  seismic 
damping  on  pure  shear  and  dilatational  modes  of  deformation.  The  trans¬ 
formation  is  to  a  set  of  generalized  coordinates  based  on  the  allowed 
element  vibrational  modes.  The  application  of  linear  transformations  on 
the  global  as  well  as  local  level  offers  the  possibility  of  further  reduc¬ 
ing  arithmetic  operation  counts,  thereby  increasing  execution  speed  of 
the  algorithm. 


(  km ) 
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Figure  1.1-4  Vertical  velocity  seismograms  for  the  basin-range-bas 
wise  linear  velocity  and  normal  surface  traction. 


2.0  GEOPHYSICAL  DATA  BASE  AND  MODELS 


Full  wave  field  seismic  calculations  are  much  more  complicated  in  3-D 
than  in  2-D.  The  reasons  are  that  3-D  geologic  data  and  2-D  surface  modeling 
and  input  processing  are  necessary,  finite  element  data  storage  and  element 
processing  time  both  increase  substantially,  and  the  amount  and  complexity 
of  output  processing  grows  manyfold.  Despite  these  added  complexities,  with 
well  designed  pre- and  post-processors  and  a  general  3-D  finite  element  solver, 
3-D  problems  can  be  reasonably  modeled,  solved  and  interpreted.  There  remains 
the  question  of  size  and  utility  for  the  scale  of  Air  Force  applications. 

An  answer  to  this  question  of  utility — whether  large-scale  3-D  simula¬ 
tions  are  practical  on  today's  supercomputers — requires  first  and  foremost  a 
good  geophysical  data  base.  In  other  words,  a  well  controlled,  practical 
geologic  model,  preferable  with  seismic  data  available  for  comparison  with 
synthetics.  Yucca  Flat,  Nevada  was  chosen  on  the  basis  of  data  collected 
from  AFOSR-sponsored  research  concerning  near-field  observations  of  signifi¬ 
cant  azimuthal  variations  in  recorded  ground  motion.  Stump  and  Reinke  (1983), 
and  far-field  observations  of  systematic  deviations  in  peak  body  wave  magni¬ 
tude  estimates  (mb) ,  Ferguson  and  Herrin  (1982) .  Other  good  sources  of  data 
are  the  detailed  geophysical  investigations  done  in  Yucca  Flat  to  insure  shot 
containment,  Cogbill  and  App  (1983). 

With  the  geologic  model  controlled  from  the  above  studies,  it  is  next  a 
matter  of  gridding  the  model,  i.e.,  coding  the  3-D  finite  element  discreti¬ 
zation,  and  running  it  on  a  CRAY.  This  is  accomplished  in  our  work  using  a 
general  purpose,  explicit  finite  element  code  called  FLEX,  consisting  of 
standard  FORTRAN  subroutines  designed  for  optimized  execution  on  scalar  as 
well  as  vector  machines.  For  the  2-D  problems  described  above,  a  standard 
CRAY-1  (one  million  words  of  memory)  was  sufficient.  However,  for  3-D  much 
more  memory  is  required  and  at  least  a  CRAY-1S  (three  million  words)  is 
necessary. 

2.1  YUCCA  FLAT  GEOPHYSICAL  MODELS 

Our  desire  to  perform  calculations  in  a  well  controlled  model  led  to 
selecting  the  north-central  part  of  Yucca  Flat,  in  the  immediate  neighborhood 
of  a  recent  underground  nuclear  explosion,  COALORA.  Extensive  near-field 
ground  motion  data  have  been  collected  from  this  shot  in  a  cooperative  effort 


(Air  Force  Weapons  Lab,  Los  Alamos  National  Lab  and  University  of  California, 
Berkeley)  to  extend  the  observational  data  base  from  contained  nuclear  explo¬ 
sions,  Stump  and  Reinke  (1983).  To  constrain  their  COALORA  propagation  model. 
Stump  and  Reinke  used  a  coarser  3-D  basin  model  developed  by  Ferguson  and 
Herrin  (1982)  from  a  variety  of  geophysical  data.  This  was  refined  in  the 
neighborhood  of  COALORA  using  geological  cross  sections  and  acoustic  logs 
provided  by  Allen  Cogbill  and  Fred  App  of  the  Containment  Group  at  LANL. 

The  resulting  propagation  model  developed  by  Stump  and  Reinke,  Fig.  2.1-1, 
is  assumed  to  be  uniformly  layered  east  and  west  of  the  fault.  Major  features 
of  the  model  include  the  Yucca  fault  dipping  approximately  72  degrees  to  the 
east  in  a  350-500  meter  layer  of  alluvium  overlying  600  meters  of  interbedded 
tuffs — all  overlying  Paleozoic  rock  at  depths  near  1  km.  Data  from  two  bore¬ 
holes  indicate  that  the  interbedded  tuffs  have  vitrified  layers  with  imped¬ 
ance  jumps  of  a  factor  of  two  or  more.  These  are  modeled  by  a  single  cap 
layer  of  Rainier  Mesa  tuff  in  Fig.  2.1-1.  The  water  table  is  placed  at  a 
depth  of  600  meters. 

Actually,  the  Rainier  Mesa  formation  dips  15-18  degrees  to  the  west  and 
is  shallower  to  the  north.  Also,  the  alluvium-tuff  interface  is  highly  irreg¬ 
ular  over  much  of  the  basin,  due  in  part  to  tectonic  movement,  Fernald  et  al. 
(1967).  In  the  COALORA  model,  smooth  depth  changes  and  3-D  variations  were 
ignored  with  respect  to  the  abrupt  change  at  the  Yucca  fault,  allowing  the 
authors  to  calculate  full-wave  synthetic  seismograms  in  the  uniformly  layered 
eastern  and  western  half space  idealizations. 

Location  of  the  shot  with  respect  to  the  Yucca  fault  is  drawn  in 
Fig.  2.1-2  and  indicates  a  smooth,  right-lateral  offset  in  the  fault  trace 
(identified  in  the  field  by  a  low  scarp  over  most  of  its  25  km  length)  within 
.5  km  southeast  of  the  shot  epicenter.  This  is  a  significant  3-D  feature, 
being  so  close  to  the  shot,  and  provides  the  central  structural  element  of 
our  3-D  geologic  model.  The  secondary  element  is  the  major  discontinuity 
between  tuff  and  Paleozoic  rock.  A  coarse  model  of  the  contact  surface  was 
provided  by  Ferguson  (1983),  based  on  studies  of  gravity,  borehole  and  seismic 
reflection  data.  Limited  resolution  of  the  data  smooths  out  any  local 
details,  including  the  Yucca  fault  which  is  known  to  have  a  vertical  offset 
of  100-200  meters  at  the  Paleozoic  contact. 


2.2  THE  2-D  COALORA  MODEL 

Before  examining  Che  3-D  model  in  Yucca  Flat,  a  2-D  plane  strain  model 
of  Che  COALORA  sice  was  investigated  to  compare  results  with  layered  half space 
calculations  in  Stump  and  Reinke  (1983),  as  well  as  recorded  seismograms. 

The  shot  was  located  in  alluvium  at  a  depth  of  280  meters,  just  above  the 
Rainier  Mesa  tuff  layer  on  the  western  side  of  Yucca  fault.  Lateral  scale 
of  the  model  covers  the  range  of  most  fielded  seismic  instruments — three 
kilometers  east  and  west  of  the  shot.  Model  depth  is  l.S  km,  400-500  meters 
into  the  underlying  Paleozoic  rock.  To  mimic  an  infinite  half space,  normal 
impedance  boundary  conditions  are  applied  to  the  two  sides  and  bottom  of  the 
finite  element  grid.  This  yields  a  simple  but  effective  absorbing  boundary 
which  is  easily  vectorized — in  contrast  to  some  ocher  boundary  formulations 
attempted. 

The  number  of  elements  in  the  model  is  chosen  to  achieve  a  given  fre¬ 
quency  resolution  for  the  slowest  body  or  surface  wave.  Assuming  negligible 
excitation  of  Rayleigh  waves  due  to  source  depth,  the  slowest  wave  is  the 
shear  body  wave  in  the  alluvium,  where  Vs  *  .98  km/s.  Assuming  10  Hz  as  the 
highest  frequency  of  interest  in  this  calculation  yields  a  minimum  shear 
wavelength  of  about  98  meters.  A  rule  of  thumb  for  explicit  finite  element 
algorithms  is  that  at  least  ten  elements  are  necessary  per  wavelength  to 
support  a  wave  with  minor  dispersion.  Therefore,  an  element  dimension  of 
9  meters  was  used,  giving  167  by  666  square  elements  in  a  1.5  by  6  km  model, 
a  total  of  111,222  elements. 

The  line  source  is  a  pressurized  two-dimensional  cavity  of  elements 
through  the  COALORA  shot  point,  280  meters  below  the  surface  in  alluvium  and 
300  meters  west  of  Yucca  fault.  Output  in  the  form  of  synthetic  seismograms 
is  obtained  by  convolving  the  calculated  velocity  response  on  the  free  surface 
(from  the  unit  pressure  step)  with  a  10  Hz  wavelet.  No  instrument  response 
is  included.  The  10  Hz  velocity  wavelet  is  a  higher  frequence  signal  than 
would  be  expected  from  shots  in  alluvium,  Haskell  (1967),  but  facilitates  the 
visual  separation  of  arrivals  in  the  synthetics. 

2.3  THE  3-D  COALORA  MODEL 

The  three-dimensional  geologic  model  in  Yucca  flat  is  contained  within  a 
3-D  rectangular  box,  1.2  km  deep  and  4.6  km  on  a  side,  centered  on  the  COALORA 
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shot  point.  Principal  features  of  the  model  are  the  Yucca  fault  with  about 

800  meters  of  right  offset  in  the  fault  trace  south  of  the  shot,  and  the 

» 
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Paleozoic  contact  surface  showing  significant  depth  variation  within  a  3  km 
radius.  The  fault  trace  is  approximated  by  a  straight  north-south  line  with 
a  linearly  offset  segment,  i.e.,  a  piecewise  linear  trace,  and  the  resulting 
fault  surfaces  dip  72  degrees  to  the  east.  The  Paleozoic  contact  is  assumed 

1 

merely  to  truncate  the  alluvium  and  tuff  basin  fill  from  below;  there  is  no 
attempt  to  model  conformable  layers  over  the  contact  surface.  A  drawing  of 
the  resulting  3-D  model  is  shown  in  Fig.  2.3-1. 

i 

When  going  from  2-D  to  3-D  modeling,  if  it  is  necessary  to  maintain  the 
same  frequency  resolution,  the  number  of  elements  goes  up  tremendously.  For 
example,  to  extend  the  6  x  1.5  km  2-D  COALORA  model  (666  x  167  ■  111,222 

elements)  to  a  6  x  6  x  1.5  ka  3-D  model  with  the  same  frequency  resolution 

would  require  666  x  111,222*  74  million  elements.  This  is  over  two  orders  of 

i 
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magnitude  larger  than  the  size  currently  feasible  with  in-core  processing  on 

the  largest  CRAY  supercomputer,  the  CRAY-1S,  with  3  million  words  of  memory. 

The  data  structure  of  FLEX,  the  finite  element  code  used  for  these  calculations. 

9 

permits  the  processing  of  500,000  3-D  elastic,  small-displacement  elements  in 

■ 

3  megawords  of  memory.  System  overhead  and  other  uncontrollable  factors  reduce 

this  number  to  slightly  less  than  400,000  elements.  Assuming  an  overall  model 

aspect  ratio  (length  or  width  to  depth)  of  5,  the  400,000  available  elements 

yield  a  block  125  x  125  x  25  -  390,625. 

9 

Clearly,  the  time  and  length  scales  of  our  3-D  model  must  be  increased 

I 

substantially  over  those  of  the  2-D  model  to  fit  into  the  CRAY-1S. 

Lowering  the  maximum  shear  wave  frequency  resolution  from  10  Hz  to  2.5  Hz 

Increases  the  element  size  by  the  same  factor,  and  decreases  the  number  of 

elements  by  4**3»64.  This  necessitates  only  a  moderate  reduction  of  the  2-D 

model  dimensions.  The  final  model  is  1.2  x  4.6  x  4.6  km — discretized  into 

125  x  125  x  25  *  390,625  cubic  elements,  37  meters  on  a  side  with  a  maximum 

frequency  resolution  of  2.6  Hz  for  shear  waves  in  the  alluvium. 

I 

The  elastic  point  source  is  represented  by  a  3-D  cavity  of  elements 
enclosing  the  shot  point.  The  cavity  is  pressurized  by  a  unit  step  at  t~0  and 

three-component  velocities  are  recorded  at  points  on  four  lines  over  the  free 
surface.  The  cavity  consists  of  three  elements  crossing  on  three  orthogonal 
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lines.  This  provides  an  approximation  intermediate  between  the  circularized 
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source  of  Che  2-0  calculation  and  che  three  dipole  forces  equivalent  to  an 
explosive  point  source,  e.g..  Love  (1944).  Lines  are  oriented  north-south, 
east-west  and  along  the  two  diagonals,  northwest-southeast  and  northeast- 
southwest.  The  source  cavity  and  output  lines  are  illustrated  in  Fig.  2.3-2. 

Synthetic  seismograms  are  calculated  by  convolving  the  recorded  veloci¬ 
ties  with  a  3.3  Hz  wavelet,  which  is  closer  to  the  expected  frequency  content 
of  the  COALORA  shot  than  the  10  Hz  wavelet  applied  in  the  2-D  calculation. 
This  frequency  is  1  Hz  higher  than  the  allowable  shear  wave  frequency  in  the 
alluvium  but,  because  P-waves  are  the  principal  phases  there  and  are  well 
modeled,  the  higher  frequency  is  used  to  enhance  the  temporal  separation  of 
these  phases  in  the  synthetics. 
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Figure  2.1-2.  Plan  view  of  Che  COALORA  site  showing  Che  Yucca 

fault,  the  6  km  east-west  line  for  the  2-D  model, 
and  the  2.3  x  2.3  km  square  for  the  3-D  model. 

For  reference,  the  line  examined  in  detail  by 
Ferguson  and  Herrin  (1982)  is  also  included. 
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Figure  2.3-1.  Yucca  Flat  3-D  model  showing  Yucca  fault  and  COALORA 
shot  location.  Model  is  4.6  km  on  a  side  and  1.2  km 
deep,  with  125  x  125  x  25  *  390,625  elements. 
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Output  lines  over  free  surface 


Source  Cavity 


2.3-2.  Illustrations  of  output  lines  over  source  cavity 
in  Yucca  Flat  model,  and  the  seven  element  3-D 
source  cavity. 


3.0  THE  FINITE  ELEMENT  CALCULATIONS 

The  calculations  described  here  were  performed  on  a  CRAY-1S  at  the  Los 
Alamos  National  Laboratory.  Computer  time  was  generously  provided  by  the  Con¬ 
tainment  Group  at  the  Lab.  Because  the  CRAY-1S  is  on  a  restricted  partition, 
it  was  not  accessible  from  our  remote  site.  Allen  Cogbill  of  the  Containment 
Group  provided  the  essential  interface  by  submitting  run  files  and  making  out¬ 
put  files  accessible.  Raw  output  data  were  transmitted  over  the  phone  between 
Los  Alamos  and  Menlo  Park  and  all  post-processing  was  done  on  an  inhouse  mini¬ 
computer. 

3.1  2-D  COALORA  CALCULATION 

The  2-D  model  and  synthetic  velocity  seismograms  are  shown  in  Figs.  3.1-1 
and  3.1-2  for  vertical  and  radial  components,  respectively.  The  seismograms 
(particle  velocity  versus  time)  are  plotted  vertically  over  each  receiver 
point  on  the  model  surface.  Amplitude  in  each  record  is  scaled  to  the  maxi¬ 
mum  first  arrival  over  the  model.  Scale  magnification  factors  are  indicated 
above  the  seismograms.  The  wavelet  has  a  duration  of  about  .18  seconds  and  a 
wavelength  of  approximately  270  meters  in  the  alluvium.  The  receiver  point 
near  ground  zero  is  offset  to  the  east  slightly,  explaining  the  radial  motion 
just  above  the  shot  in  Fig.  3.1-2. 

Within  a  range  of  1  km,  the  first  arrival  is  seen  to  be  the  direct  wave 
through  the  alluvium.  Beyond  1  km,  a  much  weaker  wave,  apparently  a  head  wave 
from  the  faster  vitrified  tuff  layer  (Rainier  Mesa  formation),  precedes  the 
direct  wave.  Just  above  the  shot,  the  vertical  record  (Fig.  3.1-1)  shows  a 
phase  shifted  wavelet — from  an  initial  zero  phase  input — caused  by  interfer¬ 
ence  between  the  direct  wave  and  that  reflected  from  and  reverberating  in  the 
tuff  layer  immediately  below  the  source.  No  significant  reflection  from  the 
1  km  Paleozoic  interface  is  evident  in  the  records  near  the  calculated  .69 
second  arrival  time.  Radial  records  indicate  a  strong  reflection  within  a 
range  of  300  meters,  arriving  at  .45  to  .5  seconds  and  comparable  in  ampli¬ 
tude  to  the  direct  wave.  This  appears  to  be  a  reflection  from  the  water  table 
at  a  depth  of  500  meters.  Comparing  vertical  and  radial  records  overall, 
radial  motion  is  stronger  and  more  coherent  than  vertical  motion,  and  records 
to  the  west  are  somewhat  more  complicated  than  those  to  the  east. 

Referring  now  to  the  radial  synthetics,  Fig.  3.1-2,  the  Yucca  fault  does 
not  appear  to  exert  a  strong  local  influence  on  the  model  response.  No  direct 
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diffractions,  attributable  to  the  fault  discontinuity  in  the  Rainier  Mesa  tuff 
j  layer  are  evident  in  the  records.  Although  certainly  present,  they  are 

obscured  by  the  direct  and  reflected  arrivals.  Differences  between  east  and 
west  records  are  probably  due  to  the  respective  depths  of  alluvium  rather  than 
to  diffractions  and  scattering  from  the  fault  itself.  Principal  differences 
|  are:  direct  waves  in  alluvium  decay  faster  to  the  west  (see  the  unmagnified 

set  of  radial  synthetics) ,  and  also  to  the  west  there  appears  a  phase  between 
1  and  1.8  km  range,  arriving  from  1.3  to  1.7  seconds  which  is  absent  to  the  east. 

Regarding  phase  velocities,  measurement  from  the  records  gives  the  direct 

) 

I  wave  in  alluvium  at  1.33  km/s,  but  the  head  waves  arriving  beyond  1  km  yield 

3.75  km/s.  The  alluvium  phase  velocity  is  2  percent  greater  than  the  wave 
speed  and  within  modeling  error  bounds.  However,  the  head  wave  velocities 
are  clearly  too  high  to  be  associated  with  the  3.25  km/s  Rainier  Mesa  tuff, 
and  are  instead  due  to  refractions  from  the  saturated  tuff  below  the  water 
table.  True  head  waves  from  the  Rainer  Mesa  tuff  are  just  noticeable  as 
first  arrivals  between  1  to  1.6  km. 

Some  limited  timing  comparisons  east  and  west  of  the  fault  can  be  made 
from  the  records  to  later  compare  with  results  from  Stump  and  Reinke  (1983). 

At  a  range  of  567  meters,  travel  time  is  435  ms  and  there  is  no  difference 
between  east  and  west — as  would  be  expected  at  this  close  range  with  no 
intervening  homogeneities  in  the  alluvium.  At  a  range  of  2.62  km,  the  timing 
discrepancy  is  56  ms,  retarded  to  the  east  and  due  to  the  depth  difference. 

The  records  are  shown  in  Fig.  3.1-3. 

3.2  3-D  COALORA  CALCULATION 

Results  of  the  3-D  calculations  are  presented  in  synthetic  velocity  seis¬ 
mograms,  Figs.  3.2-1  through  3.2-3,  for  vertical,  radial  and  transverse  compo¬ 
nents.  The  synthetics  consist  of  42  individual  seismograms  equally  spaced  on 
lines  passing  over  the  shot.  The  west-east  and  south-north  (global)  lines  are 
4.5  km  long  (111  meter  spacing),  and  the  diagonal  lines  are  6.36  km  (157  meter 
spacing),  oriented  southwest-northeast  and  northwest-southeast.  Three  sets  of 
synthetics  are  displayed  for  each  line,  with  magnifications  of  1,  20  and  2000. 

The  high  magnification  is  used  to  examine  the  weak  refracted  phases.  Duration 
of  the  synthetics  is  1.98  seconds  after  pressurization  of  the  cavity. 
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3.2.1  Vertical  Velocity  Synthetics 


Vertical  velocity  synthetics  are  shown  in  Fig.  3.2-1  for  the  four  lines. 
Because  the  predominant  input  frequency  is  a  third  of  that  used  in  the  2-D 
synthetics,  wavelets  are  three  times  longer,  hence  less  distinct  and  interfere 
to  a  greater  degree.  The  strongest  motion  is  found  over  the  shot  point  at 
around  .4  seconds,  caused  by  constructive  interference  of  the  direct  wave  and 
reflected  wave  from  the  layer  of  Rainier  Mesa  tuff.  Referring  to  the  unmagni¬ 
fied  synthetics,  ground  motion  decays  rapidly  with  range — due  in  part  to 
scaling  by  peak  interference  amplitude  rather  than  peak  direct  wave  amplitude. 

At  a  magnification  of  20  the  synthetics  clearly  show  the  direct  wave  as 
first  arrival  out  to  .9  km,  at  a  measured  phase  velocity  of  1.52-1.58  km/s  on 
all  lines.  Beyond  a  range  of  1.35  km  the  first  arrival  is  a  refraction  and 
the  phase  velocity  apparently  depends  on  the  depth  of  alluvium  under  the  line 
observed.  On  the  S-N  line,  phase  velocity,  V  is  3.1  km/s  to  the  south  and 
3.26  km/s  to  the  north.  On  the  W-E  line,  V  is  3.1  km/s  to  the  east  and  3.26 
km/s  to  the  west.  From  Fig.  2.3-1,  the  south  and  east  ends  lie  on  . 5  km  of 

alluvium  on  the  downthrown  side  of  the  fault,  and  the  north  and  west  ends  are 

over  .35  km  of  alluvium  on  the  upthrown  side.  Therefore,  calculated  phase 
velocities  are  lower  on  the  down  side  (east  and  south)  and  higher  on  the  up 
side  (west  and  north).  By  overlaying  the  synthetics  so  that  east  is  plotted 
over  south  and  west  over  north,  we  find  that  all  arrivals  are  coincident, 
independent  of  range. 

On  the  diagonal  lines,  also  in  Fig.  3.2-1,  the  refracted  phases  are  less 
distinct.  Approximate  measurements  of  slope  on  the  synthetics  at  a  magnifi¬ 
cation  of  20  give  phase  velocities  of  the  first  discernable  arrivals  as  3.25- 

3.32  km/s  on  the  northwest  and  southwest  ends,  respectively.  Consequently, 
phase  velocities  on  the  diagonals  do  not  show  the  trend  nor  the  disparity 

noted  above  between  up  and  down  sides  of  the  fault.  On  the  southwest  end  of 

the  SE-NE  line  the  phase  velocity  is  seen  to  increase  significantly  between 
-2.5  and  -3.2  km.  This  is  attributed  to  a  refracted  wave  in  the  shallow 
Paleozoic  rock  underlying  the  southwest  corner  of  the  model.  Fig.  2.3-1. 

3.2.2  Head  Waves 

To  see  the  earliest  arriving  signals  in  the  finite  element  model,  we 
must  examine  the  synthetics  under  very  high  magnification.  At  2000X  there  is 
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a  clear  first  arrival  which  precedes  all  of  the  much  stronger  arrivals  des¬ 
cribed  above.  Referring  to  the  south-north  line  in  Fig.  3.2-1,  phase  veloc¬ 
ities  are  found  to  be  3.0  km/s  on  the  south  end  (down  side)  and  2.7  km/s  on 
the  north  end  (up  side) — in  contrast  to  the  trend  observed  above  for  later  but 
stronger  arrivals.  Physically,  these  first  arrivals  should  be  head  waves  from 
the  Rainier  Mesa  tuff  layer. 

Calculating  the  theoretical  head  wave  arrival  time  at  a  range  of  1.8  km 
gives  .76  seconds  at  the  north  end  and  .83  seconds  at  the  south  end.  Measured 
times  from  the  seismograms  at  this  range  are  .74  seconds  and  .83  seconds  north 
and  south,  respectively ,  within  2-3  percent  of  the  theoretical.  Therefore,  in  terms 
of  arrival  times  the  phase  is  a  refraction  in  the  tuff  layer,  but  phase  veloc¬ 
ities  are  8  percent  and  16  percent  below  the  theoretical  value  of  3.25  km/s. 
Discrepancies  in  V  are  caused  in  part  by  the  3-D  finite  element  discretization, 
which  only  uses  three  elements  through  the  thickness  to  represent  the  high 
speed  tuff  layer.  This  is  clearly  adequate  to  model  arrival  times  within  a 
few  percent,  but  measured  phase  velocities,  which  depend  on  the  angle  of  inci¬ 
dence  of  the  wave  from  below,  are  more  sensitive  to  the  coarse  discretization 
in  the  layer. 


3.2.3  Precursor  Waves  from  Grid  Dispersion 

Further  examination  of  the  2000X  synthetics  shows  that  just  above  the 
shot,  the  earliest  arrival  time  is  .045  seconds,  compared  to  the  theoretical 
body  wave  time  of  .15  seconds.  This  discrepancy  is  associated  with  the  finite 
element  discretization  in  general,  and  the  3-D  source  region  in  particular. 

When  the  source  cavity  is  pressurized,  in  addition  to  the  radial  modes 
of  an  ideal  spherical  cavity,  extraneous  deformational  modes  are  generated. 
Most  of  these  propagate  at  the  appropriate  body  wave  speed  and  decay  rapidly 
with  range — however,  weak  precursors  caused  by  grid  dispersion  will  always 
propagate  one  element  per  timestep.  For  a  timestep  of  .0066  seconds,  with 
seven  elements  between  the  free  surface  and  the  shallowest  pressurized  ele¬ 
ment,  the  theoretical  arrival  time  is  .046  seconds  (versus  .045  seconds  mea¬ 
sured  from  the  seismograms)  for  an  apparent  body  wave  speed  of  5.64  km/s. 

The  precursors  may  therefore  be  identified  in  the  2000X  seismograms  as 
the  first  phase  observed  immediately  above  the  source.  Included  in  these 
precursors  is  the  so-called  hourglass  mode  which  is  suppressed  using  a  spe¬ 
cial  hourglass  stiffness  parameter  in  the  finite  element  algorithm.  This 


spurious  node  is  only  significant  (under  high  magnification)  in  the  immedi¬ 
ate  neighborhood  of  the  source.  It  will  always  be  present  to  some  degree  and 
is  most  prevalent  when  only  a  few  elements  are  contained  in  the  source  cavity, 
one  pressurized  element  being  the  worst  case.  Hourglassing  will  be  discussed 
further  in  the  theoretical  development. 

3.2.4  Radial  Velocity  Synthetics 

Radial  velocity  components  on  the  global  west-east  and  south-north  lines 
are  shown  in  Fig.  3.2-2.  The  unmagnified  synthetics  on  global  lines 
indicate  that  radial  motion  is  more  coherent  and  decays  less  rapidly  with  range 
than  vertical  motion  on  the  same  line.  Comparing  the  magnified  radial  records 
to  the  verticals  discussed  in  3.2-1  shows  essentially  the  same  first  arrival 
times  and  phase  velocities  over  each  end  of  the  lines.  Within  .9  km  of  the 
shot,  radials  on  the  diagonal  lines  exhibit  a  stronger,  slightly  earlier  direct 
arrival.  The  implication  is  that  vertical  and  radial  motions  are  well  corre¬ 
lated,  i.e.,  surface  ground  motions  arise  from  the  same  incident  body  wave 
phases. 


3.2.5  Transverse  Velocity  Synthetics 

Transverse  synthetics  on  the  global  lines  in  Fig.  3.2-3  show  relatively 
weak  arrivals  with  peak  motions  25-30  percent  of  the  radial  component.  The 
transverse  arrivals  appear  to  correlate  with  vertical  and  radial  components 
of  motion.  Careful  comparison  shows  that  transverse  peaks  are  retarded  by 
.15  seconds  on  both  ends  of  the  W-E  line,  while  on  the  S-N  line  peaks  are 
coincident  on  the  south  end  and  retarded  by  .15  seconds  on  the  north  end. 
Motion  asymmetry  on  the  S-N  line  is  clearly  visible  in  the  20X  synthetics, 
with  earlier,  stronger  transverse  velocities  on  the  south  end.  The  asymmetry 
might  be  attributed  to  wave  interaction  at  the  fault  crossing,  south  of  ths 
shot — probably  diffractions  from  the  obliquely  oriented  discontinuity  in  the 
Rainier  Mesa  tuff  layer.  No  such  discontinuity  is  crossed  on  the  north  end 
and  the  transverse  arrivals  are  much  weaker  there.  On  the  W-E  line  the  arriv¬ 
als  are  weak  on  both  ends,  although  the  east  line  crosses  the  fault  perpen¬ 
dicular  to  the  tuff  discontinuity  and  the  west  line  does  not  cross  at  all. 
Therefore,  waves  incident  on  the  fault  trace  at  an  oblique  angle  appear  to 


excite  stronger  transverse  motions  by  diffraction  from  the  edge  of  tuff  layer 
discontinuities  across  the  fault. 

Referring  now  to  transverse  arrivals  on  the  diagonal  lines,  also  in 
Fig.  3.2-3,  there  is  a  marked  difference  in  arrival  times  between  synthetics 
on  the  diagonals  and  those  on  the  global  lines.  At  IX  the  transverse  velocity 
is  clearly  higher  on  the  southeast  and  northeast  ends  of  the  lines.  These 
ends  both  cross  the  fault  trace  obliquely  and  appear  to  support  the  hypothesis 
that  diffractions  from  oblique  crossings  of  the  discontinuity  are  responsible 
for  the  stronger  transverse  arrivals. 

The  transverse  phases  on  the  diagonals  also  show  a  large  time  delay  in 
comparison  to  phases  on  the  globals.  If  diffraction  is  the  principal  source 
of  these  phases, then  the  retarded  arrivals  should  agree  with  theoretical 
timing  from  the  source  cavity  to  the  fault  edge  and  up  to  the  surface.  The 
first  arrival  at  the  fault  is  a  weak  head  wave  through  the  Rainier  Mesa  tuff 
layer  at  .12  seconds,  and  travel  time  of  the  resulting  diffracted  wave  to  the 
free  surface  above  the  shot  is  .37  seconds— giving  a  total  time  delay  of  .49 
seconds.  Measured  delay  in  the  20X  synthetics  is  .45  seconds  for  the  weak 
precursors  and  about  .5  seconds  for  the  true  arrival.  This  strongly  suggests 
that  diffractions  from  the  fault  are  the  principal  source  of  transverse  phases 
on  the  diagonal  lines,  which  is  further  supported  by  the  skewed  appearance  of 
the  synthetics,  i.e.,  the  apparent  source  is  shifted  to  the  right,  consistent 
with  the  direction  and  approximate  range  of  the  diffracting  fault  edge  from 
the  source  cavity. 

The  strong  asymmetry  in  arrival  times  on  either  end  of  the  diagonal  lines 
also  appears  to  be  caused  by  the  fault.  Measured  phase  velocities  on  the  east 
(down)  side  of  the  fault  (right  side  of  the  synthetics)  are  1.55-1.7  km/s  on 
both  lines,  while  on  the  west  (up)  side  they  are  3. 3-3. 5  km/s.  Therefore 
arrivals  to  the  east  propagate  through  the  slow  alluvium,  and  to  the  west  they 
propagate  through  the  faster  tuffs.  The  reason  that  no  wave  energy  is  seen  in 
the  tuffs  to  the  east  is  that  only  weak  secondary  diffractions  (diffractions 
of  diffractions)  from  the  up  side  of  the  fault  edge  can  couple  into  the  down- 
thrown  layers.  In  contrast,  diffracted  energy  from  the  discontinuity  is  cou¬ 
pled  directly  into  the  layer  on  the  west  side.  This  is  illustrated  in 
Fig.  3.2-4. 

Comparing  transverse  motion  on  the  diagonals  to  that  on  the  globals  shows 
earlier  and  stronger  arrivals  on  the  global  W-E  and  S-N  lines.  The  discrepancy 
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is  caused  in  pare  by  a  one-half  element  offset  of  the  global  lines  from  the 
true  shot  center,  causing  an  erroneous  transverse  component  which  is  strongest 
near  the  source.  The  diagonal  lines  pass  directly  over  the  shot  center. 
Quantitatively*  the  diagonal  line  transverse  motions  are  20-30  percent  of  the 
global  line  motions  and  the  half  element  offset  does  not  seem  to  account  for 
the  full  discrepancy. 
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Figure  3.1-1.  The  2-D  COALORA  model  and  vertical  velocity  seismograms 
calculated  on  the  free  surface.  The  source  cavity  is 
pressurized  at  t  *  0.  The  static  water  level  is  indi¬ 
cated  by  SWL. 
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Figure  3-1-3  Vertical  velocity  seismograms  for  a  10  Hz  wavelet  at  600  m  and  2619  m  ranges 
east  and  west  of  Yucca  Fault  showing  calculated  time  delays. 
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Figure  3.2-1 . 


Vertical  velocity  synthetic  seismograms  with  magnifications 
indicated  on  right.  Peak  velocity  is  .367  x  10“?  m/s  and 
the  upper  set  is  for  diagonal  lines  while  the  lower  is  for 
global  lines. 
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4.0  THEORETICAL  ANALYSIS  OF  THE  EXPLICIT  WAVE  SOLVER 


Large-scale  simulations  of  time  domain  seismic  waves  in  heterogeneous 
3-D  geology  require  discrete  numerical  models  with  hundreds  of  thousands  to 
millions  of  degrees  of  freedom.  Finite  difference  or  finite  element  methods 
reduce  these  models  to  a  large  system  of  ordinary  differential  equations  in 
time.  The  simplest  explicit  integration  scheme  (lumped  mass  approximation) 
provides  a  natural,  as  well  as  practical  and  efficient  algorithm  for  stepping 
the  system  forward  in  time.  The  scheme's  utility  follows  from  the  property 
that  all  degrees  of  freedom  are  uncoupled  during  each  time-step.  This  explicit 
decoupling  is  a  direct  expression  of  the  wave-like,  i.e.,  hyperbolic  nature 
of  the  continuum  problem,  namely  that  disturbances  at  one  point  take  a  finite 
amount  of  time  to  influence  another  point.  In  the  discrete  model  the  decou¬ 
pling  is  explicitly  implemented  by  choosing  a  time-step  less  than  the  shortest 
transit  time  between  any  two  nodes.  The  clear  advantage  over  implicit  algo¬ 
rithms  is  that  no  large  coupled  system  of  algebraic  equations  is  necessary. 

The  resulting  computer  code  loops  are  readily  vectorized  for  optimal  process¬ 
ing  on  pipelined  machines. 

Explicit  integration  methods  are  not  without  certain  disadvantages,  how¬ 
ever.  One  is  the  inherent  difficulty  in  modeling  viscous  material  damping  in 
a  manner  consistent  with  the  damping  behavior  of  a  physical  medium.  Geologic 

media  are  generally  characterized  by  the  quality  factors,  Q  and  Qa  for 

ot  p 

dilatational  and  shear  dominated  motion,  respectively.  However,  it  is  not 
clear  using  the  conventional  numerical  formulation  how  to  solve  for  the  model 
damping  coefficients  as  functions  of  Q  and  Qa  purely  on  the  basis  of  nodal 
field  quantities.  In  other  words,  given  Q's  over  the  model,  how  can  viscous 
damping  coefficients  for  the  orthogonal  components  of  nodal  velocity  be 
determined?  One  approach  suggested  by  structural  models  is  to  assume  so- 
called  Rayleigh  damping,  e.g..  Bathe  (1978)  where  the  model  damping  coefficient 
is  a  linear  combination  of  the  coefficients  multiplying  displacement  and 
acceleration  vectors,  i.e.,  the  stiffness  and  mass  matrices.  In  its  general 
form  of  implementation  this  is  a  poor  approximation  for  a  continuum.  However, 
in  terms  of  a  more  natural  set  of  local  basis  vectors  (rather  than  the  nodal 
basis),  useful  alternatives  to  Rayleigh  damping  may  be  derived. 

Another  disadvantage  to  explicit  methods  is  the  presence  of  a  class  of 
unrestrained  deformational  modes,  so-called  hourglass ing  modes.  These  are 
by-products  of  the  approximate  spatial  integration  scheme  used  in  calculating 


nodal  field  quantities,  i.e.,  internal  restoring  forces.  The  modes  are 
associated  with  linear  variations  in  the  field  between  discrete  nodes  and  can 
be  selectively  and  artificially  restrained  after  isolating  them  from  the  local 
grid  deformations,  e.g.,  Flanagan  and  Belytschko  (1981).  Hourglass  modes  are 
usually  unimportant  components  of  the  wave  field  except  near  abrupt  changes 
in  boundary  or  loading  conditions,  or  discontinuous  2-D  and  3-D  geometry  like 
corners  and  cavities.  In  the  physical  model  these  same  discontinuities  are 
strong  sources  of  diffracted  waves,  implying  a  correspondence  between  diffrac¬ 
tions  and  hourglassing.  The  degree  to  which  hourglass  suppression  influences 
calculated  diffracted  fields  is  unknown. 

In  this  theoretical  section  of  the  report  a  discrete  analyses  will  be 
performed  to  investigate  aspects  of  explicit  wave  solution  theory  relevant  to 
the  above  comments.  These  include  Che  vector  formulation  of  an  explicit 
finite  element  algorithm  in  terms  of  a  very  useful  modal  element  basis.  This 
basis  extends  the  theory  of  hourglass  suppression  to  selective  damping  of 
dilatational  and  shear  modes  and  is  derived  by  diagonalizing  the  element  mass 
and  stiffness  matrices. 

4.1  THE  FINITE  ELEMENT  EQUATIONS— NODAL  BASIS 

To  optimize  the  processing  of  discrete  numerical  models  on  pipelined 
machines,  it  is  necessary  to  fully  vectorize  the  algorithm  while  minimizing  the 
number  of  arithmetic  operations.  In  the  following,  a  straightforward  set  of 
vector  algebraic  equations  is  derived  to  accomplish  this,  based  on  a  somewhat 
different  view  of  the  explicit  finite  element  method.  The  approach  is  to 
reduce  the  governing  element  equations  to  their  canonical  form  by  transforma¬ 
tions  to  element  deformation  basis  vectors  which  diagonalize  the  mass  and 
stiffness  matrices. 

The  finite  element  grid  used  to  discretize  the  model  is  assumed  to  be 
Cartesian,  i.e.,  composed  of  uniform  rectangular  elements.  Global  Cartesian 
grids  provide  the  simplest  formulation  possible  and  eliminate  the  arithmetic 
for  intermediate  element  transformations  mapping  irregular  hexahedra  to 
Cartesian  boxes.  Of  course  this  also  eliminates  the  ability  to  conform 
smoothly  to  irregular  (non-Cartesian)  boundaries.  However,  one  virtue  of 
large-scale  models  is  that  the  grid  is  sufficiently  dense  to  allow  good  resolu¬ 
tion  with  only  stepwise  approximations  to  irregularities.  In  2-D  this  is 
analogous  to  the  use  of  many  discrete  picture  elements  (pixels)  to  represent 
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a  graphics  Image.  If  che  model  can  resolve  a  wave  pulse  of  length ,  A,  then  it 
can  also  resolve  a  material  inhomogeneity  (object)  with  similar  characteris¬ 
tic  length.  As  a  consequence,  either  the  highest  frequency  or  smallest 
object  to  be  represented  determines  a  consistent  grid  spacing. 

The  finite  element  model  consists  of  a  three-dimensional,  rectangular 
domain  partitioned  by  a  uniform  Cartesian  grid  into  identical  cubic  elements. 
Fig.  4.1-1.  The  unknown  displacements  are  to  be  determined  at  each  grid 
point  or  node.  Displacements  within  elements  are  given  by  a  shape  function 
which  interpolates  nodal  values  and  maintains  continuity  between  elements. 

The  simplest  shape  function  for  a  Cartesian  element  is  a  3-D  version  of 
the  trapezoidal  rule.  It  can  be  derived  formally  from  the  local  3-D  Taylor 
series  expansion  by  neglecting  curvature  terms,  i.e.,  terms  with  second  or 
higher  derivatives  in  one  of  the  space  variables.  This  function  gives 

u(x°  +  ?,  t)  =  uC  +  C,  Uc  +  G-  uC  4-  r  uC  +  uC 

-  ~  ~  1  ^  2  ~x2  3  ~x3  1  2  x^x^ 

+  51?3  ~X]x3  +  ?2?3  -x2x3  +  ^ 1C2^3  -x1x,x3  (1) 

T 

where  u  ■  (u^,u2,u3)  is  the  displacement  vector,  x  and  Q  are  global  and 
local  coordinate  vectors,  respectively,  and  c  superscripts  refer  to  element 
centroid  values.  The  24  components  of  the  superscripted  vector  on  the  right 
can  be  solved  in  terms  of  the  24  nodal  displacements  of  the  element,  giving 

u(x°+5,  t)  *  D(C)  6(xC,  t)  (2) 

T 

where  6  *  (u^, . . .  *u^g,  u^, . . .  ,u2g,  u3^,...,u3g^  is  the  nodal  displace¬ 
ment  vector  and  subscripts  refer  to  the  coordinate  and  node  number, 
respectively. 

With  an  analytical  expression  available  for  displacements  within  the  ele¬ 
ment,  the  theory  of  elasticity  gives  strains  and  stress  as 
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where  C  is  the  matrix  of  constitutive  constants  and  B  is  the  strain- 
displacement  matrix,  found  by  differentiating  (2)  with  respect  to  4^, 
and  substituting  into  the  strain  expression, 


The  final  step  applies  to  principle  of  virtual  work  to  calculate 
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generalized  elastic  and  inertial  forces,  f  and  f  corresponding  to  the  nodal 
displacements,  6.  For  a  virtual  displacement,  6  the  elastic  work  is 

/gTOdv  -  /5TBTCBSdv 
v~  ~  v~ 

and  the  inertial  work  due  to  constant  mass  density,  p  is 

/u^(pu)dv  *  /5^D^pD6dv 
v~ 

/•T  E  »T  I 

External  work  of  the  elastic  and  Inertial  nodal  forces  are  6  f  and  -6  f  , 
respectively.  Equating  these  to  the  integrals  and  solving  gives 

fE  -  K<S,  K  =  /BTCBdv  (5) 

v 

f1  -  -M<$,  M  =  p/DTDdv  (6) 

v 

where  K  and  M  are  the  symmetric  stiffness  and  mass  matrices  and  the  integrals 

are  over  element  volume.  Therefore,  the  disc  '.te  system  of  equations  governing 
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motion  of  the  single  finite  element  is,  from  aynamic  equilibrium,  f  *  f  +  g, 

hence  M6  *  -K6  +  g  (7) 

where  g  represents  the  influence  of  adjacent  elements,  external  forces  and  damp¬ 
ing.  The  24  x  24  matrices  are  easily  evaluated  and  K  is  found  to  be  fully 
populated  while  M  consists  of  three  8x8  blocks  on  the  main  diagonal. 

4.2  THE  FINITE  ELEMENT  EQUATIONS— MODAL  BASIS 

Element  stiffness  and  inertia  properties  can  be  defined  with  respect  to 
any  complete  set  of  displacement  modes.  The  nodal  basis  used  above  is  ideal 
for  deriving  the  system  but  not  necessarily  for  optimal  computational 
efficiency.  For  the  latter  purpose,  it  is  necessary  to  carefully  examine  the 
structure  and  symmetries  of  the  system  of  24  equations  governing  a  single 


element.  A  formal  diagonalization  of  the  mass  and  stiffness  matrices  would 

yield  the  forma,  _ 

M  ■  (8) 

K  -  QXQT  (9) 

where  P  and  Q  are  matrices  of  eigenvectors  and-«  and  K  are  diagonal  eigenvalue 
matrices.  This  is  an  effective  approach  but  the  required  eigenvalue-eigenvector 
analysis  is  non-trivial.  Fortunately,  a  simpler  alternative  exists. 

The  alternative  is  to  diagonalize  the  matrices  by  inspection,  on  the 
basis  of  displacement  modes  of  a  unit  cube.  It  is  natural  to  start  with  the 
mass  matrix,  M,  because  of  its  block  structure.  The  three  8x8  blocks  are 
each  associated  with  one  component  of  nodal  displacement  so  only  one  block 
need  be  considered,  the  others  following  by  a  permutation  of  indices.  To 
proceed,  recall  that  formal  diagonalization  of  square  matrices  and  modal 
analysis  of  discrete  structures  share  the  same  underlying  mathematics,  differ¬ 
ing  only  in  viewpoint.  Therefore,  diagonalization  of  the  block  matrices  may 
be  achieved  by  a  single  modal  analysis  of  the  cube  under  the  constraint  that 
only  one  component  of  displacement  is  involved.  These  modes  consist  of  a 
rigid  body  translation  and  seven  deformations.  The  modal  shapes  are  repre¬ 
sented  by  their  eight  ordered  nodal  displacements  in  vector  form  and  all 
modes  are  orthogonal  (i.e.,  scalar  products  vanish). 

The  seven  deformation  modes  are  of  two  kinds.  The  first  kind  includes 
the  three  modes  familiar  from  classical  elasticity  theory — one  extension  (or 
compression)  and  two  shears  in  say,  the  5^.  direction.  The  second  kind  con¬ 
sists  of  four  so-called  hourglass  modes.  These  are  less  familiar  because 
they  represent  bending  deformations  which  are  higher  order  and  therefore 
vanish  in  the  usual  mathematical  limiting  process  with  respect  to  modes  of 
the  first  kind.  Modes  of  the  second  kind  are  encountered  in  the  couple- 
stress  theory  of  elasticity,  Sternberg  (1968)  which  provides  a  natural 
theoretical  basis  for  their  further  study. 

The  complete  set  of  mass  matrix  modes  for  the  cube  is  plotted  in 
Fig.  4.2-1,  where  each  is  represented  by  its  proportions1  "odal  displacement 
vector,  S*,  n  ■  1-8  and  i  ■  1  for  the  direction.  The  rigid  body  mode  is 

a  vector  with  all  elements  equal.  Modes  of  the  first  kind  are  given  by 
vectors  with  elements  proportional  to  the  ?2 *  ant*  £3  n°dal  coordinates. 
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respectively,  in  the  cube-centered  Cartesian  coordinate  system,  Fig.  4.1-1. 
Modes  of  the  second  kind  may  be  derived  by  an  orthogonalization  process 
starting  from  the  four  known  orthogonal  mode  vectors  (1  rigid  body  and  3 
first  kind),  e. g. ,  Gram-Schmidt  orthogonalization,  or  by  simple  experimenta¬ 
tion.  They  are  also  available  in  the  literature,  e.g.,  Flanagan  and  Belytschko 
(1981). 


To  diagonalize  the  mass  matrix  the  block  modal  or  eigenvector  matrix  is 
formed  as  P  *  (S^...Sg)  from  which  the  similarity  transformation,  P  in  (8) 
becomes 


(10) 


The  P^/2 3  and  P  matrices  are  unitary  (orthogonal)  meaning  that  P-^  *  P1.  The 


diagonalized  mass  matrix  can  therefore  be  written  as 


T 

P  MP 


Diagonalization  of  the  stiffness  matrix  proceeds  by  first  applying  the 
similarity  transformation  for  the  mass  matrix,  P  as 

K  *  PTKP 

Detailed  analysis  and  results  from  the  interim  report,  Wojcik,  et  al.  (1982) 
show  that  K  is  considerably  simpler  than  K,  but  not  yet  diagonal.  The  number 
of  nonzero  elements  is  reduced  from  576  in  K  to  39  in  K  because  of  the  greatly 
simplified  stress  state  for  the  element  modes  in  Fig.  4.2-1.  The  18  off- 
diagonal  terms  arise  from  coupling  of  the  and  components  of  the 

modes  of  the  first  kind,  i.e.,  coupling  between  the  three  extensional  modes 
and  between  the  six  shear  modes.  Modes  of  the  second  kind  are  uncoupled. 
Therefore  it  is  necessary  to  diagonalize  the  reduced  system  as 

*:  -  rtkr 

where  matrix  R  is  determined  by  generalizing  the  2-D  results  obtained  in  the 
interim  report.  In  particular,  the  9  modes  of  the  first  kind  (3  modes  in  3 
directions)  can  be  decomposed  into  3  Poisson  stretch  modes,  3  pure  shear 
modes  and  3  rigid  rotation  modes.  These  stiffness  modes  are  drawn  in 
Fig.  4.2-2  for  the  direction  and  can  be  written  in  terms  of  linear  combi¬ 
nations  of  mass  modes  of  the  first  kind,  resulting  in  an  expression  for  R 
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consisting  of  constants  in  the  nonzero  locations  of  K.  Therefore,  the  diago¬ 
nal  stiffness  matrix  becomes 

T  T  T 

X  -  KP  KPR  -  Q  KQ 

In  terms  of  the  above  factorization  of  the  mass  and  stiffness  matrices, 
the  equations  for  a  simple  finite  element  become 

PAPT6  -  -PRXRTPT5  +  g  (11) 

Therefore,  to  solve  for  the  nodal  accelerations,  it  is  only  necessary  to  pre- 

-1  T 

multiply  through  by  the  inverse  of  the  factored  mass  matrix,  i.e.,  P jl  P  ,  giving 

5  -  PX^RKrVs  +  P^_1PTg  (12) 

This  is  a  new  form  of  the  consistent  mass  equation  for  nodal  acceleration — 
new  in  terms  of  the  simple  factorization  of  M~^K  multiplying  6.  In  most 
practical  applications  to  continuua,  the  consistent  mass  matrix,  M  is 
replaced  by  a  diagonal  lumped  mass  approximation,  M  to  allow  a  trivial  inver¬ 
sion  of  M,  hence  a  solution  for  6.  This  yields 

6  -  M~1PRXRTPT6  +  ft_1g  (13) 

neglecting  the  external  force  vector,  g  and  comparing  the  consistent  mass 
equation,  (12)  to  the  lumped  mass  approximation,  (13)  shows  that  in  terms  of 
arithmetic  operations  there  is  no  advantage  to  mass  lumping  in  the  above 
formulation. 

4.3  VISCOUS  DAMPING 

There  remains  the  question  of  a  consistent  viscous  damping  term  in  the 
finite  element  equations  of  motion.  To  make  the  damping  term  explicit  in  our 
formulation  the  force  vector,  g  is  replaced  by 

g-+  -  C5  +  g 

The  unknown  damping  matrix,  C  is  factored  as 

C  -  SCST 

where  C  and  S  are  its  eigenvalue  and  eigenvector  matrices.  The  intent  is  to 
damp  extenaional  and  shearing  motions  of  the  medium,  so  it  is  natural  to  con¬ 
sider  the  case  where  S  •  Q,  i.e.,  where  the  damping  matrix  has  the  same 
eigenvectors  as  the  stiffness  matrix.  The  eigenvalues  of  C  therefore  eoasist 
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of  Che  damping  coefficients  of  the  individual  modes  in  the  stiffness  bases— 
the  rigid  translations  and  rotations,  Poisson  extensions,  pure  shears  and  the 
stress-couple  or  hourglass  modes.  It  follows  that  each  of  these  modes  of  a 
finite  element  can  be  selectively  damped  by  choosing  an  appropriate  constant 
for  each  eigenvalue  of  C  in  Che  stiffness  basis. 

The  resulting  damped  finite  element  equations  become,  from  (11), 

P*PT6  +  QCQT5+  QXQT<$  =  g  U*) 

In  practice  the  only  modes  for  which  damping  data  is  available  are  the  exten- 
sional  and  shear  modes.  In  order  to  evaluate  the  damping  constants  in  (14) 
from  the  given  quality  factors,  and  Qg  (or  alternatively,  given  percentages 
of  critical  damping) ,  exact  solutions  for  the  natural  frequencies  of  the 
Poisson  extension  and  pure  shear  modes  are  necessary.  These  can  be  obtained 
from  a  vibrational  analysis  of  Che  two  deformation  modes  in  Fig.  4.2-2. 


5.0  CONCLUSIONS 


The  conclusions  drawn  from  this  research  effort  concern:  large-scale 
numerical  modeling  of  seismic  waves  on  supercomputers;  a  geophysical  data 
base  adequate  to  justify  3-D  simulations  and  interpretations  of  very  near¬ 
field  seismic  events;  and  a  theoretical  study  of  specializations  of  the 
explicit  finite  element  algorithm  used  for  the  calculations.  The  primary 
conclusion  is  that,  given  the  three-dimensional  geologic  database  available 
from  past  and  ongoing  investigations  in  Yucca  Flat,  Nevada  Test  Site,  the 
size  of  feasible  3-D  computational  models  on  the  CRAY-1S  appears  to  be 
adequate  to  simulate  full  elastic  wave  fields  and  interpret  arrivals  for  com¬ 
parison  with  existing  3-D  ground  motion  data.  The  most  significant  short¬ 
coming  is  lack  of  an  analytical  implementation  of  the  3-D  pressure  source. 
Instead,  a  relatively  crude  discretization  of  the  source  cavity  is  pressur¬ 
ized,  producing  some  artificial  asymmetry  in  the  outgoing  pressure  pulse. 

Thus  the  discrete  source  in  3-D  requires  further  investigation  and  characteri¬ 
zation.  In  two  dimensions,  much  higher  grid  density  is  available  for  modeling 
so  that  discretization  errors  are  not  significant  in  large-scale  2-D  applica¬ 
tions. 

Synthetic  seismograms  from  the  3-D  simulation  of  the  COALORA  event  at 
Yucca  Flat  indicate  that  a  significant  source  of  transverse  motion  on  radial 
lines  through  the  source  is  diffraction  from  abrupt  changes  in  geology.  In 
our  model,  diffractions  were  observed  from  the  discontinuity  in  the  Rainier 
Mesa  layer  across  the  Yucca  fault.  This  suggests  that  local  geologic  inhomo¬ 
geneities  near  the  source  will  act  as  simple  diffraction  sources — effectively 
providing  secondary  seismic  sources  in  the  free-field.  This  is  in  contrast 
to  the  mechanisms  of  free-field  scattering  and  strain  release  which  can  also 
contribute  to  transverse  motion.  Future  large-scale  calculations  of  plane 
waves  through  a  detailed  3-D  inhomogeneous  model  in  Yucca  Flat  would  enable 
us  to  quantify  the  strength  of  the  diffracted  transverse  motion. 

Regarding  diffracted  fields,  their  representation  by  explicit  finite 
element  methods  requires  further  investigation.  The  reason  is  that  the 
pressurized  source  and  material  discontinuities  produce  hourglass  modes  which 
are  artificially  suppressed  during  computations.  However,  we  observe  that 
diffraction  and  hourglassing  are  closely  associated  because  they  are  produced 


by  similar  features,  namely  discontinuities  in  loading  or  geometry.  There¬ 
fore,  what  is  the  effect  of  houglass  suppression  on  real  diffracted  waves? 

The  answer  to  this  question  must  come  from  a  theoretical  study  of  the  finite 
element  algorithm. 

Successful  time-domain  simulations  in  3-D  are  clearly  feasible  with 
pipelined  supercomputers.  Optimal  processing  still  requires  careful  tailor¬ 
ing  of  the  algorithm,  however,  to  vectorize  inner  code  loops  and  eliminate 
nonessential  arithmetic.  An  effective  approach  to  simplify  the  implementa¬ 
tion  and  theory  is  to  use  dense  Cartesian  grids,  typical  of  finite  difference 
methods,  to  discretize  the  model.  The  simplest  case  is  the  uniform  grid  con¬ 
sisting  of  cubic  elements — used  exclusively  in  the  present  study.  Theoretical 
analysis  of  the  cubic  finite  element  shows  that  the  application  of  heuristic 
eigenvalue-eigenvector  concepts,  i.e.,  a  change  of  basis,  yields  a  closed 
form  diagonalization  of  the  element  mass  and  stiffness  matrices.  This  pro¬ 
vides  a  drastic  reduction  in  the  number  of  floating  point  operations  required 
to  evolve  the  element  forward  in  time. 

In  addition  to  the  arithmetic  simplifications  provided  by  the  eigenvalue- 
eigenvector  element  modal  analysis,  the  resulting  displacement  basis  admits 
selective  damping  of  the  individual  vibrational  modes  of  the  element.  Con¬ 
sistent  viscous  damping  has  always  been  a  problem  with  explicit  finite  elements 
but  the  modal  basis  on  cubic  elements  appears  to  yield  an  effective  approach, 
worthy  of  further  investigation.  The  modal  basis  also  allows  a  consistent 
mass  formulation  with  the  same  number  of  arithmetic  operations  as  the  lumped 
mass  approximation,  as  well  as  an  exact  integration  to  produce  the  stiffness 
matrix.  For  comparison,  the  FLEX  code  used  for  the  present  calculations 
(representative  of  state-of-the-art  explicit  finite  element  theory)  is 
formulated  with  lumped  masses  and  approximate  single  point  Gaussian  integra¬ 
tion.  All  of  these  theoretical  simplifications  derive  from  the  cubic  element 
and  the  application  of  linear  algebra  techniques  to  the  element  matrices. 

Regarding  hourglass  suppression  and  its  Influence  on  diffractions,  an 
examination  of  the  element  modes  shows  12  so-called  hourglass  modes  for  a 
cube.  These  modes  are  involved  with  bending  of  the  element  due  to  moments  or 
couples  acting  on  the  element  faces.  However,  the  hourglass  modes  of  a  cube 
are  not  complete,  hence  overly  stiff,  because  curvature  terms  are  excluded 
from  the  element  shape  function.  By  including  12  curvature  terms  in  the 


Taylor  series  (2  for  each  face)  and  12  additional  generalized  displacements* 
a  consistent  bending  theory  can  be  derived.  The  appropriate  generalized  dis¬ 
placements  are  expected  to  be  the  12  edge  rotations.  Rotations  are  usually 
included  in  plate  bending  theories  using  Hermite  polynomials,  but  a  similar 
inclusion  in  3-D  continuum  theories  appears  to  be  new.  We  observe  chat 
bending  modes  are  associated  with  the  couple-stress  formulation  in  elasticity 
theory,  e.g.,  Mindlin  (1963),  Sternberg  (1968).  Couple-stress  theory  has 
been  thoroughly  studied  but  is  of  limited  utility;  however,  the  present  work 
indicates  that  it  may  provide  a  further  analytical  framework  for  assessing 
the  effect  of  hourglass  control  on  diffractions. 
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